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1  Overview  on  the  research  program 

The  objective  of  this  project  is  to  develop  and  analyze  efficient  and  robust  nu¬ 
merical  methods  for  simulation  of  mathematical  models  of  nonlinear  phenomena 
such  as  nonlinear  conservation  laws,  advection-dominated  flows,  and  other  related 
problems  where  nonlinear  effect  such  as  shocks,  fronts,  and  contact  discontinu¬ 
ities  pose  significant  difficulties  for  traditional  numerical  methods. 

The  main  thrust  of  the  research  program  for  the  period  July  2012  to  June  2015 
was  to  explore  and  investigate  the  approximation  properties  of  a  technique  that  we 
call  entropy  viscosity  method.  The  main  stabilization  mechanism  in  this  method  is 
a  nonlinear  dissipation  proportional  to  the  local  size  of  an  entropy  production.  We 
believe  that  the  entropy  viscosity  idea  is  an  important  conceptual  breakthrough. 
In  the  three  year  period  of  the  current  project  we  have  made  progresses  in  the 
following  directions: 

1.  Robust  explicit  viscous  method.  We  have  constructed  a  first-order  ap¬ 
proximation  technique  for  nonlinear  scalar  conservation  equations  that  is 
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maximum  principle  satisfying  on  any  mesh,  any  space  dimension  and  for 
any  flux.  The  approximation  is  explicit  in  time  and  done  with  continuous 
finite  elements;  the  flux  just  needs  to  be  Lipschitz  continuous. 

2.  Error  estimates.  We  have  proved  error  estimates  for  the  above  first-order 
method.  The  methods  converges  strongly  to  the  unique  entropy  solution 
with  a  convergence  rate  which  is  at  least  0(h*)  irrespective  of  the  space 
dimension. 

3.  Second-order  method.  We  have  proposed  an  accurate  and  maximum  prin¬ 
ciple  preserving  method  for  nonlinear  scalar  conservation  equations  using 
the  so-called  flux  correction  technique  (FCT).  The  method  combines  the 
first-order  technique  and  the  entropy- viscosity  methodology. 

4.  Viscous  regularization  of  the  Euler  equations.  We  have  investigated  vis¬ 
cous  regularization  techniques  for  the  compressible  Euler  equations  using 
fundamental  results  from  thermodynamics.  A  general  class  of  viscous  regu¬ 
larization  that  is  compatible  with  all  known  generalized  entropy  inequalities 
has  been  identified. 

5.  Positivity  in  the  Euler  equations.  We  have  started  to  analyze  the  positivity 
properties  of  explicit  viscous  schemes  using  continuous  finite  elements  for 
the  compressible  Euler  equations.  We  have  obtained  an  explicit  positivity¬ 
preserving  scheme  and  shown  that  positivity  does  not  depend  on  the  speed 
of  sound.  The  technique  does  not  involve  any  Riemann  solvers  (approxi¬ 
mate  or  exact). 

6.  Entropy-viscosity  LES.  We  have  investigated  the  idea  of  using  the  no¬ 
tion  of  entropy-viscosity  as  a  LES  technique  for  the  incompressible  Navier- 
Stokes  equations  at  high  Reynolds  numbers. 

7.  High  order  schemes  for  other  models:  We  are  developing  a  new  sec¬ 
ond  order  methods  for  the  so-called  Mean  Filed  Games  (MFG).  The  MFG 
model  is  a  coupled  nonlinear  system  which  emerges  in  financial  and  crowd 
dynamics  models. 


2  Entropy  viscosity 

We  review  in  this  section  some  elementary  facts  about  the  entropy  viscosity  method¬ 
ology.  Nonlinear  scalar  conservation  equations  can  all  be  put  into  the  following 
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general  form: 


dtu{x,  t)  +  V-/(w(x,  t))  =  0,  x  G  fl,  t  >  0,  (1) 

with  u\t=0  =  u0  and  appropriate  boundary  conditions.  The  scalar  initial  boundary 
value  problem  has  a  unique  entropy  solution  which  satisfies  an  additional  set  of 
differential  inequalities 

dtE{u)  +  V-F(u)  <  0,  (2) 

for  any  pairs  E(u)  and  F(u)  such  that  E  is  convex  and  F'(u)  =  E'{u)f'{u).  The 
function  E  is  called  entropy  and  F  is  the  associated  entropy  flux.  For  convex 
fluxes  (i.e.,  if  /  is  convex)  in  one  space  dimension  it  is  known  that  one  entropy 
pair,  for  example  the  one  generated  by  E(u)  =  \u2,  is  enough  to  select  the  unique 
entropy  solution.  Physical  systems  have  at  least  one  entropy  pair  and  the  entropy 
inequality  (|2j)  is  the  mathematical  form  of  the  second  principle  of  thermodynam¬ 
ics.  It  is  expected  that  the  auxiliary  inequality  (|2])  serves  as  a  selection  criteria 
and  guarantees  convergence  of  the  numerical  approximation  to  the  correct  physi¬ 
cal  solution  of  the  nonlinear  system.  Therefore,  it  is  desirable  (and  necessary)  to 
somehow  incorporate  the  entropy  dissipation  (j2j)  in  a  numerical  scheme. 

A  traditional  way  of  selecting  the  entropy  solution  of  ([I)  consists  of  adding 
viscous  dissipation 


dtue  +  V-f(ue)  -  V-(eVtte)  =  0,  (3) 

where  e  >  0  and  it  can  be  shown  in  general  that  ue  —t  u  when  e  —?  0.  The 
parameter  e  is  usually  taken  to  be  proportional  to  the  local  mesh  size  when  con¬ 
structing  numerical  approximation  of  ([3])  and  this  limits  the  convergence  rate  to 
first-order  at  most.  The  use  of  artificial  viscosity  to  solve  nonlinear  conservation 
equations  has  been  pioneered  by  Neumann  and  Richtmyer  and  popularized  later 
by  Smagorinsky  for  LES  purposes  and  by  Ladyzenskaja  for  theoretical  purposes 
in  the  analysis  of  the  Navier-Stokes  equations.  However,  the  early  versions  of  ar¬ 
tificial  viscosities  were  overly  dissipative  and  the  interest  for  development  of  such 
methods  has  faded  over  the  years.  The  research  direction  shifted  towards  Dis¬ 
continuous  Galerkin  and  Finite  Volume  methods,  where  upwinding  and  limiters 
have  been  shown  to  be  efficient  and  to  yield  high-order  accuracy.  Up  to  a  few 
exceptions  slope/flux  limiting  is  a  one-dimensional  concept  that  does  not  easily 
generalize  to  unstructured  meshes  in  higher  dimensions.  Moreover,  the  theoretical 
understanding  of  the  stability  and  convergence  of  limiters  is  currently  restricted  to 
uniform  grids  and  scalar  equations  in  one  space  dimension.  For  the  above  reasons 
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and,  among  other  things,  the  fact  that  artificial  viscosities  are  easy  to  implement, 
the  interest  for  artificial  viscosity  has  lately  been  revived  in  the  DG  literature  and 
in  the  Continuous  Galerkin  (CG)  literature  as  well. 

It  is  usually  argued  in  the  literature  that  good  artificial  viscosities  can  be  com¬ 
puted  from  measures  of  the  local  regularity  of  the  solution  or  from  local  residuals 
of  the  PDE.  We  propose  to  take  a  slightly  different  route  by  considering  the  local 
residual  of  an  entropy  equation  to  construct  an  artificial  viscosity. 

e~  ch2\dtE{u)  +  V-F(u)\.  (4) 

One  immediate  consequence  of  this  choice  is  that  the  viscosity  is  proportional  to 
the  entropy  production,  which  is  known  to  be  large  in  shocks  and  to  be  zero  in 
contact  discontinuities.  As  a  result,  this  strategy  makes  an  automatic  distinction 
between  shocks  and  contact  discontinuities,  and  this  subtle  distinction  cannot  be 
made  by  any  of  the  two  classes  of  methods  mentioned  above.  We  also  think 
that  using  the  residual  of  the  conservation  equation  may  be  less  robust  than  using 
the  entropy  residual.  This  argument  is  based  on  the  observation  that  consistency 
requires  the  residual  of  the  PDE  to  converge  to  zero  in  the  distribution  sense  as  the 
mesh-size  goes  to  zero,  whereas  the  very  nature  of  entropy  implies  that  the  entropy 
residual  converges  to  a  Dirac  measure  supported  in  the  shocks.  This  implies  that 
the  entropy  residual  focuses  far  better  in  shocks  than  the  PDE  residual,  and  it  is 
in  this  sense  that  we  claim  that  the  PDE  residual  is  less  reliable  than  the  entropy 
residual.  This  can  be  better  understood  by  considering  the  simple  case  of  the  one- 
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Figure  1:  Entropy  viscosity  in  log  scale  for  the  Burgers  equation  at  t  —  0.25 

dimensional  Burgers  equation  dtu  +  \dxu2  =  0  with  initial  data  u(x,  0)  =  1  if 
x  <  0  and  u{x,  0)  =  0  otherwise.  The  entropy  solution  is  u(x,  t)  =  1  —  H(x  — 
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where  H  the  Heaviside  function.  One  observes  that  the  residual  of  the  equation 
is  zero,  whereas  the  entropy  residual  is  an  negative  measure:  \dtu 2  +  ^du3  = 
—j^S(x  —  |f),  where  6  is  the  Dirac  measure.  This  effect  is  very  well  illustrated 
on  the  one-dimensional  Burgers  equation  over  the  interval  (0, 1)  with  initial  data 
sin(7nr),  see  Figure [T]  We  show  in  this  figure  the  entropy  viscosity  computed  by 
using  either  the  quadratic  entropy  E(u)  =  (first  and  second  panel  from  the 
left)  or  E(u)  =  u  (third  and  fourth  panels).  Choosing  E(u)  =  u  corresponds  to 
using  the  residual  of  the  PDE  for  the  viscosity.  These  tests  show  that  choosing 
E(u )  =  \u2  leads  to  better  focusing  of  the  viscosity  and  is  more  robust  with 
respect  to  the  multiplier  c  (see  definition  (Q). 

We  now  illustrate  the  entropy  viscosity  concept  on  the  compressible  Euler 
equations  by  using  the  physical  entropy  to  construct  the  entropy  viscosity.  We 
show  in  Figure  [2]  the  density  field  at  t  =  2.86  and  the  viscosity  field  at  t  =  4 
for  the  classical  wind  tunnel  problem  with  a  forward  facing  step  at  Mach  3.  We 
observe  that  the  viscosity  focuses  very  well  and  there  is  almost  no  viscosity  in  the 
contact  discontinuity  that  develops  at  the  top  of  the  flow. 


(a)  Density  at  t  =  2.86  (b)  Entropy  viscosity  flux  at  t  =  4 


Figure  2:  Wind  tunnel  with  a  step  at  Mach  3,  Pi  approximation. 

We  show  in  Figure[3]the  density  field  at  t  =  0.2  for  the  double  Mach  reflection 
problem  at  Mach  10.  We  observe  again  that  there  is  almost  no  viscosity  in  the 
contact  discontinuity  that  develops  from  the  first  triple  point;  moreover  the  jet 
develops  the  standard  instability  at  the  bottom  of  the  domain. 

We  have  also  developed  two  different  techniques  to  enforces  boundary  condi¬ 
tions  on  curved  boundaries,  which  is  not  an  easy  task  for  nonlinear  conservation 
equations.  We  are  also  using  the  entropy  viscosity  framework  to  construct  a  goal 
oriented  refinement  strategy.  The  method  is  well  developed  by  now.  We  have 
solved  various  benchmark  problems  with  finite  elements,  spectral  elements  and 
Fourier.  We  show  typical  results  in  Figure  |4j 
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(a)  Density  (b)  Entropy  viscosity 

Figure  3:  Double  Mach  reflection,  t  =  0.2,  Pi  approximation. 


Figure  4:  Supersonic  around  a  cylinder  in  a  tunnel.  Pi  approximation.  Density 
field  (left)  and  entropy  viscosity  (right). 

3  Achievements  in  the  period  July  2012- June  2015 

We  describe  in  this  section  what  we  have  achieved  during  the  period  July  2012- 
June  2015. 


3.1  A  robust  viscous  method 

It  is  traditional  to  stabilize  numerical  methods  for  solving  transport  equations  or 
scalar  conservation  equations  by  adding  artificial  dissipation.  In  the  continuous 
finite  element  literature  the  viscosity  is  added  by  augmenting  the  Galerkin  formu¬ 
lation  with  the  weak  form  of  the  operator  —V-{vhV)  where  vh  is  scalar- valued 
and  Vh  =  c/3h.  Here  c  is  an  adjustable  constant,  /?  >  0  is  proportional  to  the  local 
wave  speed  and  h  is  a  measure  of  the  local  mesh  size.  The  first  obstacle  one  runs 
into  when  using  this  approach  is  that  of  a  proper  definition  of  the  meshsize  h  on 
non-uniform  anisotropic  meshes.  Although  many  clever  and  reasonably  well  jus¬ 
tified  ideas  have  been  proposed  to  address  this  non-trivial  issue,  to  the  best  of  our 
knowledge,  none  of  them  has  yet  lead  to  a  provable  maximum  principle  holding 
for  every  nonlinear  flux  and  every  mesh,  assuming  piecewise  linear  approximation 
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of  course.  Once  a  meshsize  h  has  been  chosen,  one  then  runs  into  the  problem  of 
choosing  the  constant  to  multiply  f3h  to  form  the  viscosity.  Although  \  seems  to 
be  a  reasonable  choice  justified  by  the  one-dimensional  analysis  on  uniform  grids, 
we  do  not  know  of  any  rational  for  tuning  this  constant  in  two  and  three  space  di¬ 
mensions  on  arbitrary  grids  besides  heuristic  arguments  and  trial  and  errors  tests. 
Even-though  it  is  possible  to  establish  well-founded  heuristic  arguments  to  tune 
the  constant,  again  we  do  not  know  of  any  argument  yielding  a  provable  maximum 
principle  for  every  nonlinear  flux  and  every  mesh.  The  last  argument  that  finally 
lead  us  to  revisit  the  first-order  viscosity  theory  is  that  it  is  not  robust  with  respect 
to  the  shape  of  the  cells.  When  trying  to  reproduce  the  one-dimensional  argument 
in  arbitrary  space  dimension  with  continuous  finite  elements  one  observes  that  the 
convex  combination  argument  can  work  only  if  f  Vipi-V<pj  dx  <  0  for  all  pairs 
of  shape  functions,  ipi,  ipj,  with  common  support  of  nonzero  measure.  This  is  the 
well-known  acute  angle  condition  assumption.  This  condition  is  easy  to  violate,  in 
particular  in  3D,  and  makes  any  method  based  on  the  bilinear  form  f  u^Vu-Vv  dx 
not  robust. 

Starting  from  the  above  observations,  in  0]  we  have  revisited  the  first-order 
viscosity  theory  from  top  to  bottom  and  introduced  a  dissipation  mechanism  based 
on  a  graph  Laplacian.  In  this  theory  the  viscosity  is  defined  cell-wise  to  ensure 
that  the  solution  is  maximum  principle  preserving.  The  method  is  explicit,  works 
with  continuous  finite  elements,  does  not  involve  any  adjustable  constants,  can 
be  defined  on  arbitrary  meshes  in  any  space  dimension,  and  is  independent  of  the 
local  anisotropy  of  the  mesh. 

3.2  Error  estimates 

We  have  analyzed  the  convergence  properties  of  the  first-order  Lagrange  finite 
element  technique  introduced  in  @]  for  scalar  nonlinear  conservation  laws.  We 
have  proved  in  J9l  that  the  error  in  the  Lf{L\)- norm  is  at  most  0(h a)  under 
the  appropriate  CFL  condition  in  any  space  dimension  and  for  any  shape-regular 
mesh  family;  the  mesh  may  be  composed  of  an  arbitrary  combination  of  sim- 
plices,  prisms,  cuboids,  etc.  The  estimate  is  established  by  using  the  technique 
of  the  doubling  of  the  variables  introduced  by  Kruskov  (1970)  and  first  used  by 
Kuznecov  (1976)  to  prove  error  estimates.  This  approach  completely  bypasses 
the  compensated  compactness  argument.  To  the  best  our  knowledge,  this  is  the 
first  time  that  a  priori  error  estimates  have  been  established  for  an  explicit  method 
using  continuous  Lagrange  finite  elements  to  approximate  nonlinear  scalar  con¬ 
servation  equations  on  nonuniform  meshes.  Similar  results  have  been  established 
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by  Cockburn  and  Gremaud  (1996),  but  the  error  estimate  therein  is  0(h» )  and  the 
algorithm  is  a  shock  capturing  streamline  diffusion  method  using  implicit  time 

3 

stepping  and  an  artificial  viscosity  scaling  like  h*  in  the  shocks. 

We  have  shown  in  @  that  the  error  in  the  L1-norm  is  directly  controlled  by 
entropy  residuals,  i.e.,  the  key  to  obtain  an  error  bound  is  to  estimate  the  entropy 
production  for  every  Kruzkov  entropy.  Using  a  graph-Laplacian-based  viscous 
dissipation  was  a  major  breakthrough  in  this  enterprise.  This  result  strengthens 
our  point  of  view  that  entropy  residuals  are  important  objects. 


3.3  Convergence  of  high-order  schemes 

In  (Hi,  we  establish  the  Instability  of  an  entropy  viscosity  technique  applied  to 
nonlinear  scalar  conservation  equations.  First-  and  second-order  explicit  time¬ 
stepping  techniques  using  continuous  finite  elements  in  space  are  considered.  The 
method  is  shown  to  be  stable  independently  of  the  polynomial  degree  of  the  space 
approximation  under  the  standard  CFL  condition.  This  technique  is  extended  to 
the  Navier-Stokes  equations  and  is  implemented  as  LES  model  in  SFEMaNS. 

In  a  sequence  of  papers  [131 12],  we  have  addressed  the  long  standing  question 
about  convergence  of  well  known  high-order  central  schemes  towards  the  entropy 
solution  of  a  scalar  conservation  law.  The  goal  was  to  find  tools  which  will  allow 
us  to  analyze  such  methods  and  prove  their  convergence.  We  proved  in  f  13lfl2ll 
a  new  maximum  principle  for  central  second  order  schemes  based  on  nonlinear 
piecewise  linear  reconstruction.  The  novelty  is  that  the  maximum  principle  proof 
is  valid  for  general  k-monotone  fluxes  and  convergence  to  the  entropy  solution  is 
proven  if  the  flux  is  strictly  convex.  Moreover,  the  second  order  reconstruction 
is  not  reduced  to  first  order  in  the  regions  of  local  extrema.  It  has  been  an  open 
problem  how  to  prove  stability  and  error  estimates  for  such  schemes  which  work 
in  practice  but  their  local  reconstruction  violates  a  regular  maximum  principle.  In 
our  most  recent  paper,  together  with  the  former  student  Orhan  Mehmetoglu,  we 
were  able  to  prove  strong  convergence  of  the  above  mentioned  numerical  methods 
towards  the  unique  entropy  solution  in  the  case  of  a  conservation  law  with  just 
bounded  initial  data  lfT4ll.  This  is  the  largest  possible  class  of  initial  conditions 
considered  in  practice  and  the  results  for  first  order  schemes  were  classical  in  this 
setup.  However,  this  was  not  known  any  of  the  second  order  schemes  based  on 
minmod-type  limiters. 
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3.4  Construction  of  a  second-order  method 


Imposing  both  the  maximum  principle  and  high-order  accuracy  in  space  is  a  diffi¬ 
cult  task  for  numerical  methods  approximating  scalar  conservation  equations.  The 
Godunov  theorem  even  asserts  that  it  is  impossible  for  one-step  linear  approxima¬ 
tion  methods  to  be  at  least  second-order  accurate  in  space  and  be  monotonicity 
preserving.  Satisfying  the  maximum  principle  and  having  high-order  accuracy  on 
arbitrary  meshes  in  any  space  dimension  are  very  desirable  properties  in  all  appli¬ 
cations,  but  these  are  two  contradicting  requirements:  (i)  on  one  hand  one  needs 
to  add  artificial  viscosity  to  guarantee  the  maximum  principle;  (ii)  on  the  other 
hand  one  needs  to  stay  as  close  as  possible  to  the  Galerkin  discretization  of  the 
problem  to  be  high-order  accurate  in  space. 

In  0,  we  have  extended  the  first-order  technique  introduced  in  01  to  an  ex¬ 
plicit  second-order  maximum  principle  preserving  numerical  method  that  works 
on  arbitrary  meshes  in  any  space  dimension  with  any  Lipschitz  flux  using  con¬ 
tinuous  Lagrange  finite  elements.  The  four  key  ingredients  are  the  first-order 
technique  of  01,  a  novel  treatment  of  the 
consistent  mass  matrix  by  Q,  a  high-order 
technique  (entropy-viscosity  method  of  Q) 
and  the  Boris-Book-Zalesak  flux  correction 
technique.  The  main  characteristics  of  the 
new  method  are:  (i)  it  is  maximum  princi¬ 
ple  preserving;  (ii)  it  preserves  the  accuracy 
of  the  entropy  viscosity  method;  (iii)  the  dis¬ 
persion  errors  induced  by  mass  lumping  are 
corrected  in  the  flux  limiting  step.  As  an  il¬ 
lustration  of  the  method  we  show  the  solu¬ 
tion  of  the  two-dimensional  scalar  conserva¬ 
tion  equation  dtu  +  V  •  f(w)  =  0  with  the 
non-convex  flux  f(w)  =  (sin  u,  cos  u)  and 
initial  data  «(x,  0)  =  I,  if  \fx^  +  y2  <  1, 

«(x,  0)  =  |,  otherwise.  It  is  a  challeng¬ 
ing  test  case  for  many  high-order  numerical  schemes  because  the  solution  has 
a  two-dimensional  composite  wave  structure.  For  example,  it  is  known  that 
some  central-upwind  schemes  based  on  WEN05,  Minmod  2,  or  SuperBee  re¬ 
constructions  converge  to  non-entropic  solutions.  The  computational  domain 
[—2,  2]x[— 2.5, 1.5]  is  triangulated  using  nonuniform  meshes  and  the  solution  is 
approximated  up  to  t  =  1  using  Pi  finite  elements.  The  graph  of  the  limited  solu- 
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Figure  5:  KPP  solution 


(a)  2208  nodes  (b)  8560  nodes 


( 


i 


(c)  34268  nodes  (d)  135841  nodes 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


tion  on  four  meshes  (2208,  8560,  34268,  and  135841  Pi  nodes,  respectively)  are 
shown  in  Figure [5}  The  helicoidal  composite  wave  is  clearly  visible. 

To  better  evaluate  the  amount  of  maxi-  Figure  6:  Cross  section,  KPP 
mum  principle  violation  by  the  entropy  vis-  problem 
cosity  solution  and  visualize  the  action  of  the 
flux  correction  technique,  we  show  in  Fig- 
ure[6]a  cross  section  along  the  line  x2  =  x\  of 
the  graph  of  the  entropy  viscosity  and  limited 
solutions  at  t  =  1  on  four  different  meshes 
(2208,  8560,  34268,  and  135841  Pi  nodes,  respectively).  The  cross  section  for 
the  entropy  solution  is  shown  in  the  left  panel  and  that  for  the  limited  solu¬ 
tion  is  shown  in  the  right  panel.  Two  regions  are  magnified  to  emphasize  the 
slight  violations  of  the  local  maximum  principle  by  the  entropy  viscosity  so¬ 
lution.  Using  the  curvilinear  abscissa  s  =  \[2x\,  the  first  region  is  the  rect¬ 
angle  ( s,u )  G  [1.2  ±  0.3]  x[|  ±  5xl0-2],  the  second  region  is  the  rectangle 
(s,u)  G  [4.7  ±  0.6]  x  [^  ±  lxlO-2].  Observe  that  in  both  cases  the  over-  and 
under-shoots  of  the  entropy  solution  are  barely  noticeable  and  seem  to  diminish 
as  the  meshes  are  refined.  Note  again  that  the  limited  solution  satisfies  the  max¬ 
imum  principle,  and  the  convergence  of  the  limited  solution  to  the  exact  solution 
seems  to  be  monotone. 

3.5  Viscous  regularization  for  Euler 

In  our  original  implementation  of  the  entropy  viscosity  method  0],  we  used  the 
Navier-Stokes  system  as  a  viscous  regularization  of  the  Euler  equations  and  we 
built  a  high-order  method  from  there.  Although  the  numerical  results  looked  good 
on  standard  benchmark  tests  (see  Figure  [2|  Figure  [3]),  extensive  numerical  exper¬ 
iments  revealed  that  the  Navier-Stokes  system  is  not  a  robust  regularization  in 
general.  Two  key  reasons  for  that  are:  (i)  the  Navier-Stokes  regularization  has  no 
mechanism  that  can  keep  the  density  positive;  (ii)  the  Navier-Stokes  regulariza¬ 
tion  is  known  to  violate  the  minimum  entropy  principle  if  the  thermal  diffusivity 
is  set  to  zero.  Therefore,  the  viscous  model  has  oscillations  that  are  not  present 
in  the  inviscid  limit.  These  observations  have  lead  us  to  investigate  this  prob¬ 
lem  carefully.  In  iflOll  we  have  identified  all  viscous  regularizations  of  the  Euler 
system  that  are  compatible  with  the  minimum  entropy  principle  and  with  some 
generalized  entropy  inequalities.  Moreover,  we  have  identified  one  single  family 
of  regularizations  that  preserves  the  positivity  of  the  density,  satisfies  a  minimum 
entropy  principle,  and  is  compatible  with  all  the  generalized  entropies  inequalities 
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as  introduced  in  Harten  et  al.  (1998). 


3.6  Positivity  for  the  compressible  Euler  equation 

In  QO  we  have  started  to  apply  the  graph  Laplacian  strategy  to  the  Euler  equations 
and  we  have  been  able  to  construct  an  explicit  first-order  method  with  continuous 
finite  elements  that  is  capable  of  preserving  the  positivity  of  the  density  and  the 
internal  energy.  Quite  surprisingly,  the  theorem  that  we  have  established  shows 
that  the  amount  of  dissipation  that  is  needed  for  this  property  to  hold  does  not 
invoke  the  speed  of  sound;  the  viscosity  is  directly  proportional  to  the  material 
velocity  only.  The  result  holds  for  every  reasonable  equation  of  state.  This  result 
goes  against  a  few  received  ideas  in  the  field.  This  work  just  started  during  Sum¬ 
mer  2014,  but  it  promises  to  be  very  fruitful  and  has  already  open  doors  on  barely 
explored  territories. 

3.7  LES  and  entropy-viscosity 

Simulating  turbulence  is  computationally  demanding  since  large  gradients  and 
eddy-phenomena  exist  in  general  at  scales  that  cannot  be  correctly  represented  by 
the  mesh.  The  small  scales  of  turbulent  flows  produce  still  smaller  ones  through 
the  coupling  of  wave  modes  via  the  action  of  the  nonlinear  term;  this  induces 
an  accumulation  of  energy  at  the  grid  scale.  The  consequence  of  this  cascade 
phenomenon  is  that  turbulent  flows  can  be  considered  as  having  singular  behaviors 
at  the  mesh  scale.  A  possible  remedy  to  this  phenomenon  that  we  proposed  in 
Q  EH  consists  of  monitoring  the  local  kinetic  energy  balance  and  introducing  a 
localized  dissipation  in  these  regions  that  is  proportional  to  the  violation  of  this 
balance.  The  deviation  from  the  local  energy  balance,  which  we  call  the  entropy 
residual, 

Z4(x,  t)  :=  dt{\n2h)  +  V-((§u£  +  ph)uh)  -  ^A(§u£)  +  ^-Vuh:Vuh  -  f-uh 

can  be  thought  of  as  an  indicator  for  local  entropy  production  in  analogy 
with  entropy  production  for  scalar  conservation  laws.  (Here  (u h,Ph)  are 
the  approximate  velocity  and  pressure  and  h  denotes  the  grid  scale.)  In  a 
resolved  flow  Dh(pc,t)  should  be  on  the  order  of  the  consistency  error  of 
the  method;  a  small  value.  In  |[8|  we  construct  a  viscosity  proportional 
to  \Dh(x,t)\.  This  so-called  entropy-viscosity  is  defined  by  ^(x, t)  := 

min  ^cmaxfi(x)||u/l(x,  t)\\,  cgfi2(x )  j  •  The  momentum  equation  is  then 
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modified  by  adding  the  term  — V-(zze(x,  i)Vu).  In  the  definition  of  z/%x.  t),  the 
constants  cmax  and  cE  are  tunable  parameters  which  depend  only  on  the  numerical 
method  and  the  geometry  of  the  mesh.  This  ensures  that  the  viscosity  never  ex¬ 
ceeds  the  first-order  upwind  viscosity.  When  the  local  grid  size  is  small  enough  so 
that  all  the  scales  are  resolved,  the  residual  \Dh(x.,  t)  |  and  zzg(x,  t  )  are  of  the  order 
of  the  consistency  error  and  the  consistency  error  times  fl(x)2,  respectively.  The 
entropy-viscosity  is  therefore  consistent,  i.e.,  vanishes  when  all  of  the  features  of 
the  flow  are  properly  resolved  at  the  grid  scale. 

We  have  done  series  of  validation  tests  on  the  above  model  in  the  pe¬ 
riodic  cube  using  a  Fourier  code  developed  at  LANL.  These  tests  are  re¬ 
ported  in  Q  To  evaluate  the  consistency  of  the  method  we  have  tested 
an  inviscid  flow  with  the  following  two-dimensional  initial  data:  u  = 
cos(87nr)  sin(87n/),  v  =  —  sin(87rx)  cos(87n/),  w  =  0.  By  construction, 
the  flow  should  remain  two-dimensional  and  laminar  and  the  total  kinetic  en¬ 
ergy  should  be  constant  for  some  time,  until  the  accumulation  of  numerical 
round-off  errors  trigger  three- 
dimensional  instabilities.  The  Eu¬ 
ler  solution  is  computed  up  to 
t  =  4  on  a  323  grid  and  up  to 
t  =  6  on  a  5123  grid  using  the 
entropy-viscosity  technique  and 
the  Smagorinsky  model  in  the 
Fourier  DNS  code.  Insets  (a)  and 
(b)  of  Figure  [7]  show  the  time  evolution  of  the  kinetic  energy  for  the  entropy- 
viscosity  solution  and  the  Smagorinsky  solution.  It  is  striking  that  the  Smagorin¬ 
sky  solution  loses  energy  fast  even  though  the  flow  is  laminar,  whereas  the 
entropy-viscosity  solution  tracks  the  DNS  solution  rather  closely.  The  conclusion 
of  these  tests  is  that  the  entropy-viscosity  model  outperforms  the  Smagorinsky 
model  since  it  is  consistent  whereas  the  Smagorinsky  model  is  not. 

Figure [8] shows  a  table  displaying  the  relative  kinetic  energy  loss  for  the  DNS, 
entropy-viscosity,  and  the  Smagorinsky  solutions  at  time  t  =  4  on  four  grids  323, 
643,  1283,  2563.  The  DNS  does  not  lose 
any  energy  at  all  resolutions  and  the  entropy- 
viscosity  solution  does  not  lose  any  significant 
amount  of  energy,  even  at  low  resolution.  The 
Smagorinsky  solution  on  the  other  hand  has 
lost  68%  of  the  energy  by  time  t  =  4  on  the 
32s  grid  and  3.3%  on  the  2563  grid.  This  test 
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confirms  that  contrary  to  the  Smagorinsky  method,  the  entropy-viscosity  method 
does  not  dissipate  energy  in  the  laminar  regions  of  the  flow.  In  Figure  [9}  we 
show  the  energy  spectra  of  simulation  runs  at  Re  ~  6500  for  various  resolu¬ 
tions.  All  runs  are  compared  against  a  resolved  DNS  run  at  resolution  2563  (called 
“No  Model”).  Each  under-resolved  simulation  is  done  using  the  entropy-viscosity 
model  and  the  Smagorinsky  model.  Note  that  the  unregularized  DNS  fail  to  cap¬ 
ture  the  correct  spectra  as  expected,  while  the  flows  regularized  with  entropy- 
viscosity  perform  significantly  better.  Note  also  that  the  entropy  viscosity  model 
is  always  closer  to  the  DNS  spectrum  than  the  Smagorinsky  model. 


Figure  9:  Energy  spectra.  Entropy-Viscosity  (EV)  vs.  Smagorinsky  (Smag). 

The  above  work  has  been  done  before  we  developed  the  graph  Laplacian  vis¬ 
cosity  theory  in  0  and  got  rid  of  cniax,  ||u/,.(x,  t)||  and  //(x)  in  the  definition  of 
the  first-order  viscosity.  An  obvious  task  ahead  of  us  is  to  extend  the  scalar  theory 
to  the  incompressible  Navier-Stokes  equations  and  develop  a  LES  model  that  is 
parameter  free. 

3.8  Lagrangian  hydrodynamics 

Through  collaborations  with  colleagues  at  Lawrence  Livermore  National  Lab¬ 
oratory,  we  have  extended  the  entropy  viscosity  methodology  to  Lagrangian 
hydrodynamics,  O-  Our  experiments  have  Figure  10:  3D  Noh  problem, 
shown  that  the  method  extends  naturally  to 
this  setting  without  any  particular  difficulty, 
thereby  proving  again  that  the  methodology 
that  we  propose  is  very  flexible.  We  show 
here  some  computations  done  by  our  student 
Vladimir  Tomov  who  has  been  be  partially 
supported  by  the  grant  (V.  Tomov  has  currently  (a)  Density  (b)  64  MPI  tasks 
postdoctoral  position  at  Lawrence  Livermore 

National  Laboratory).  Figure  [TO]  shows  the  3D  Noh  implosion  test.  The  field 
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shown  is  the  left  panel  is  density.  The  mesh  subdivision  among  the  64  parallel 
tasks  that  were  used  for  this  computation  is  shown  in  the  right  panel.  The  method 
is  stabilized  by  using  the  entropy  viscosity  technique. 

3.9  High-order  schemes  for  other  models 

Together  with  Vladimir  Tomov  we  are  developing  a  new  second  order  scheme 
for  approximating  nonlinear  systems  called  Mean  Field  Games  (MFG).  This  type 
of  models  have  been  emerging  recently  in  financial  applications  and  modeling 
crowd  dynamics.  The  MFG  system  couples  a  forward  Hamilton- Jacobi  equation 
and  a  backward  transport  equation.  Currently,  there  are  now  known  second  order 
numerical  approximations  for  such  models.  Together  with  Tomov,  we  have  been 
working  on  the  design  of  a  second  order  approximation  of  the  inviscid  MFG  model 
and  we  are  preparing  a  paper  on  that  lfT5l. 


3.10  Mentoring 


The  work  described  above  has  been  done  in  collaboration  with  various  gradu¬ 
ate  students  and  post-doctoral  researchers.  More  specifically  the  work  done  on 
the  graph-Laplacian  viscosity  has  been  done  with  the  visiting  assistant  profes¬ 
sor  Murtazo  Nazarov.  The  contribution  of  M.  Nazarov  to  the  work  described  in 


§3.1|-§3.4|was  important.  M.  Nazarov  visited  TAMU  between  201 1  and  2014.  He 
now  holds  a  tenure-track  position  in  Uppsala,  Sweden.  The  work  described  in 
§3.4|was  also  done  with  the  help  of  the  graduate  student  Yong  Yang  (expected  to 
graduate  in  2016).  The  work  on  LES  was  done  in  collaboration  with  Adam  Lar- 
ios  (visiting  assistant  professor  (2011-2014))  and  former  graduate  student  Travis 
Thompson.  Adam  Larios  has  now  a  tenure-track  position  at  the  university  of  Ne¬ 
braska.  Travis  Thompson  (graduated  in  2013)  has  been  research  assistant  at  the 
university  of  Tennessee  Joint  Institute  for  Computational  Sciences  at  Oak  Ridge 
National  Laboratory;  he  is  now  research  assistant  at  Rice  University.  Additional 
work  on  the  use  the  entropy  viscosity  in  Lagrangian  hydrodynamics  has  been  done 
by  the  former  graduate  student  Vladimir  Tomov.  He  graduated  in  2014  and  now 
is  research  assistant  at  Lawrence  Livermore  National  Laboratory.  Some  of  the 
funds  from  the  201 1/2015  AFOSR  grant  were  used  to  support  these  post-doctoral 
researchers  and  graduate  students. 


14 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


3.11  Personnel  Supported  During  Duration  of  Grant 

1 .  Jean-Luc  Guermond,  Professor,  Texas  A&M  University. 

2.  Bojan  Popov,  Professor,  Texas  A&M  University. 

3.  Murtazo  Nazarov,  Visiting  Assistant  Professor,  Texas  A&M  University. 
Now  tenure-track  at  Uppsalla  Univ. 

4.  Vladimir  Tomov,  Graduate  student,  Texas  A&M  University.  Now  post-doc 
at  LLNL. 

5.  Yong  Yang,  Graduate  student,  Texas  A&M  University. 

6.  Manuel  Quezada  ,  Graduate  student,  Texas  A&M  University. 
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